############################################################################
# RE main function
# (simplified for only the specification required in the bootstrap)
############################################################################

RE_func_sim = function(data,data_ext=0,quiet=TRUE,min0=NA,max0=NA,min0_mg=NA,max0_mg=NA,Nsim=4000) {

  dt = copy(data)
  dt_ext = copy(data_ext)
  
  # Prepare data for simulated densities and CDF
  min_val=min(data$coef,min0,na.rm=TRUE)
  max_val=max(data$coef,max0,na.rm=TRUE)
  dt_sim = data.table(id=seq(1,Nsim),x=seq(min_val,max_val,length.out=Nsim))
  summary(dt_sim)
  
  ############################################################################
  # Prior 
  ############################################################################
  
  dt[,ones:=1]
  
  # Prior Mean
  w_name="w"
  ww = dt[,.SD,.SDcols=w_name]
  setnames(ww,"w")
  ols_mean = lm(coef ~ z ,data=dt,weights=ww$w)
  dt[,mu0 := predict(ols_mean)]
  
  # Prior Variance
  dt[,beta_tilde:=coef-mu0]
  dt[,yn:=beta_tilde^2-se^2]
  dt[,se2:=se^2]
  dt[,var0:=wtd.mean(yn,weights=w)]
  if (mean(dt$var0)<=0) {
    dt[,var0:=wtd.var(coef,weights=w)-wtd.mean(se2,weights=w)*shrink]
  }
  
  dt=dt[,.SD,.SDcols=c("id","z","sd_z","coef","se2","mu0","var0","w")]

  ############################################################################
  # Signal 
  ############################################################################
  
  var_noise = wtd.mean(dt$se2,weights=dt$w)
  var_noise
  var_FE = wtd.var(dt$coef,weights=dt$w)
  var_FE
  mean_FE = wtd.mean(dt$coef,weights=dt$w)
  mean_FE
  
  ############################################################################
  # Posterior distribution
  ############################################################################
  
  # Posterior means
  dt[,rho:=var0/(var0+se2)]
  dt[,mu1:=mu0*(1-rho) + rho*coef]
  summary(dt)
  
  # Posterior variances
  dt[,var1:=var0*(1-rho)]
  summary(dt)
  
  ############################################################################
  # Marginal effects
  ############################################################################
  
  s_mg=NA

  #### Paste the mu1 and var1 of coefficients to data_ext 
  dt_ext = merge(dt_ext,dt[,.(id,coef,mu0,var0,mu1,var1)],by="id",all=TRUE)
  dt_ext[, mg:=coef]
  summary(dt_ext)
  
  # #### Generate data for simulated distribution
  min_mg=min(dt_ext$mg,dt_ext$mu0,min0_mg,na.rm=TRUE)
  max_mg=max(dt_ext$mg,dt_ext$mu0,max0_mg,na.rm=TRUE)
  
  #### Characterize posterior distribution
  dt_ext[,ones:=1]
  q = c(0.1,0.25,0.5,0.75,0.9)
  cdf1_mg= corrected_cdf(data=dt_ext,w_name="w",Nsim=Nsim+1,min_val=min_mg, max_val=max_mg,quiet=quiet)
  quan1 =corrected_quantiles(mat_cdf=cdf1_mg,q=q)
  v =c(quan1,corrected_mean(dt_ext,w_name="w"),corrected_var(dt_ext,w_name="w"),nrow(dt_ext))
  s_mg = as.matrix(v,ncol=1)
  rownames(s_mg)=c(paste0("Percentile ",c(10,25,50,75,90)),"Mean","Variance","N")
  
  ############################################################################
  # Return
  ############################################################################
  
  return(list(dt=dt,dt_ext=dt_ext,table_mg=s_mg))
  
}

